Cell-type-specific transcription factor analysis during HSC development

1.1 Data introduction

The single-cell RNA-sequencing data including 6 sequential cell populations during HSC development (aortic endothelial cells (AECs) and HSC-primed hemogenic-endothelial cells (HECs) from E10 aorta-gonad-mesonephros (AGM) region, T1 and T2 pre-HSCs from E11 AGM, E12 & E14 fetal liver(FL) HSCs, and HSCs from adult bone marrow).

The TF activity estimation followed with the workflow at metaRegulon vignette. Here, we don’t repeat the GRN inference and TF activity estimation steps, but showing the following steps such as cell-type-specificity analysis.

1.2 Loading data and generating objects.

The raw expression matrix and TF activity data can be downloaded from Zenodo.

suppressMessages(suppressWarnings(library(PUIC)))
suppressMessages(suppressWarnings(library(scATFR)))
suppressMessages(suppressWarnings(library(metaRegulon)))
atfr <- readRDS("D:/HSC_Autophagy/Markdown/5.metaRegulon/metaRegulon/inst/extdata/HSC_development/atfr.rds")

Next, we can extract the TF activity object act_sce that deposit in alternative assays from the atfr object for the following analysis.

act_sce <- altExp(x = atfr, e = "viper")
class(act_sce)

## [1] "SingleCellExperiment"
## attr(,"package")
## [1] "SingleCellExperiment"

Now, lets take a look at the values in the act_sce object:

assay(act_sce)[1:5,1:5]

##               E10.0_AEC_1_1 E10.0_AEC_1_2 E10.0_AEC_1_3 E10.0_AEC_1_4
## 2700081O15Rik     2.6691509     3.4384611      3.519365     2.8035412
## Adnp              0.9074099     1.5317953      1.168628     0.5021225
## Adnp2             0.3717352     0.0147019      1.855357     2.8388917
## Aebp1            -0.4412220    -1.2217474     -0.107998    -0.6754781
## Aebp2             9.1200368    10.7359366     10.183974    10.5986185
##               E10.0_AEC_1_5
## 2700081O15Rik      3.189519
## Adnp               1.525927
## Adnp2              2.124061
## Aebp1             -1.506020
## Aebp2              9.507078

We can see that the values in act_sce have negative values that may affect the following analysis steps, here we performed a liner normalization step to scales the values into 0 to 10.

assay(act_sce) <- t(apply(assay(act_sce),1,function(x){10*(x-min(x))/(max(x)-min(x))}))
assay(act_sce)[1:5,1:5]

##               E10.0_AEC_1_1 E10.0_AEC_1_2 E10.0_AEC_1_3 E10.0_AEC_1_4
## 2700081O15Rik      3.256182     4.5867276      4.726653      3.488614
## Adnp               5.732479     7.3491813      6.408844      4.683081
## Adnp2              1.642584     0.7662129      5.284269      7.698447
## Aebp1              8.256700     6.4896524      9.011094      7.726363
## Aebp2              4.952973     6.2213735      5.788111      6.113586
##               E10.0_AEC_1_5
## 2700081O15Rik      4.156174
## Adnp               7.333986
## Adnp2              5.943830
## Aebp1              5.846082
## Aebp2              5.256781

Then, we can also generate another object exp_sce including the expression of TFs obtained form atfr object.

exp_sce <- atfr[row.names(act_sce),]
exp_sce

## class: SingleCellExperiment 
## dim: 831 154 
## metadata(2): scATFR atfr_version
## assays(2): counts logcounts
## rownames(831): 2700081O15Rik Adnp ... Zxdc Zzz3
## rowData names(0):
## colnames(154): E10.0_AEC_1_1 E10.0_AEC_1_2 ... Adult_HSC_3739
##   Adult_HSC_3740
## colData names(3): sample stage sizeFactor
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(1): viper

2.1 TF expression level

library(scater)
exp_sce <- runPCA(x=exp_sce, exprs_values="logcounts") 
exp_sce <- runTSNE(x = exp_sce, exprs_values="logcounts")
plotReducedDim(exp_sce, dimred="TSNE",colour_by = "stage")

We noticed that the order of stages is according to alphabet, but this is not our desire order. To custom the order of the stages, we can use reLevel() function:

exp_sce <- reLevel(x = exp_sce, colable="stage", level_order=c("AEC", "HEC", "T1_pre_HSC", "T2_pre_HSC", "FL_HSC", "Adult_HSC"))
p <- plotReducedDim(exp_sce, dimred="TSNE",colour_by = "stage")
p

The parameter colable indicate the column label used for order adjustment, while the level_order indicate the updated order.

2.2 TF activity level

Next, we use TF activity to perform dimensional reduction analysis

act_sce <- reLevel(x = act_sce, colable="stage", level_order=c("AEC", "HEC", "T1_pre_HSC", "T2_pre_HSC", "FL_HSC", "Adult_HSC"))
act_sce <- runPCA(x=act_sce, exprs_values="atfr") 
act_sce <- runTSNE(x = act_sce, exprs_values="atfr")
q <- plotReducedDim(act_sce, dimred="TSNE",colour_by = "stage")
q

To compare the dimensional reduction between TF expression and TF activity and save the results, one can use the following codes:

ggplot2::ggsave(filename ="tf_expression_vs_activity_tSNE_res.pdf" ,plot = p+q, device = "pdf", 
       path = "./", width = 150, height =50, dpi = 600, limitsize = FALSE, units ="mm")

2.3 Indicated TF expression vs activity

Here we introduce the demonstration of TFs at expression or activity level with tSNE plots. Since the dimensional reduction results of TF expression and TF activity is quit different, we use the tSNE result of TF activity as the coordinate. Firstly, we add logcounts matrix from exp_sce into object act_sce:

assay(act_sce,"logcounts") <- assay(exp_sce,"logcounts")[row.names(act_sce),]

Then, we can compare the TF expression and TF activity of the key TFs as follow:

i <- "Hes1"
p1 <- plotReducedDim(object = act_sce,dimred = "TSNE", colour_by = i, by_exprs_values = "logcounts", point_size=0.5, point_alpha=1)
p2 <- plotReducedDim(object = act_sce,dimred = "TSNE", colour_by = i, by_exprs_values = "atfr", point_size=0.5, point_alpha=1)
p1+p2

3.1 Differential activity analysis between two adjacent stages

Here we performed the differential analysis between AECs and HECs at TF activity level using R package scran:

library(scran)
AEC_HEC_diff_act <- findMarkers(x = act_sce, groups=colData(act_sce)$stage, restrict=c("AEC","HEC"))
AEC_HEC_diff_act

## List of length 2
## names(2): AEC HEC

The results of differential activity analysis is a list of Data Frames.

HEC_pos_diff_act <- as.data.frame(AEC_HEC_diff_act$HEC)
HEC_pos_diff_act$Direction <- ifelse(HEC_pos_diff_act$FDR<0.05,ifelse(HEC_pos_diff_act$summary.logFC>0,'up','down'),'not_sig')
HEC_pos_diff_act$Direction <- factor(x = HEC_pos_diff_act$Direction,level=c('up','down','not_sig'))
head(HEC_pos_diff_act)

Now, we can use volcano plot to visualize the results as follow:

library(ggplot2)

vol <- ggplot(data = HEC_pos_diff_act, mapping = aes(x = summary.logFC, y=-log10(p.value), color=Direction))+
  geom_point(shape=16)+
  scale_color_manual(values = c("red","blue","gray"))+
  theme_bw()
vol

3.2 cell-type-specific TF comparison between expression and activity level

Here we use findAllCTSRegulons() to calculate the cell-type-specific TFs for all 6 stages at regulon level and expression level.

cts_act <- findAllCTSRegulons(x = act_sce, clust_name = 'stage', use_altExp = FALSE, use_assay = "atfr", nperm = 10000, ncores=4)
cts_exp <- findAllCTSRegulons(x = exp_sce, clust_name = 'stage', use_altExp = FALSE, use_assay = "logcounts",nperm = 10000, ncores=4)

Please note that if you want to calculate the cell-type-specific TFs for only few stages, you can use findCellTypeSignature() to calculate them. For example, To calculate the TFs that specifically activated in AEC:

AEC_cts_act <- findCellTypeSignature(x = act_sce, pos_clust="AEC", clust_name = 'stage',use_altExp = FALSE,nperm=10000,ncores=4)

Next, we want to identify TFs that only showed cell-type-specificity at activity level but not expression level:

library(purrr)
library(dplyr)
cts_act_list <- split(cts_act,cts_act$stage)
cts_exp_list <- split(cts_exp,cts_exp$stage)
cts_act_only <- purrr::map2(.x = cts_act_list, .y = cts_exp_list,.f = function(x,y){
  x %>% filter(!(regulon %in% y$regulon))
})
cts_act_only <- do.call(rbind,cts_act_only)
head(cts_act_only)

3.3 Visualization of cell-type-specific TFs

To explore the distribution these cell-type-specific TFs, we utilize the RadViz figure to visualize the distribution of these signature genes. Firstly, we calculate the average activity of each TF in every stage using function summarizeClusterValue().

mean_act_score <- summarizeClusterValue(object = act_sce, clust_name = "stage", gene_name = unique(cts_act_only$regulon), use_assay = "atfr")
mean_act_score <- data.frame(gene_name = row.names(mean_act_score), mean_act_score)
head(mean_act_score)

Then, we use the hexagonal density plot to visualize the distribution of cell-type-specific TFs as follow:

cts_act_frm <- left_join(x = cts_act_only, y = mean_act_score,by=c("regulon"="gene_name"))
stages <- c("AEC", "HEC", "T1_pre_HSC", "T2_pre_HSC","FL_HSC","Adult_HSC")
cts_act_frm$stage <- factor(cts_act_frm$stage,levels=stages)
rv1 <- RadvizEnrichPlot(x = cts_act_frm, color_by='stage', anchors = stages, optim_anchor = F, recenter = stages, text_size = 4, plot_type = "hexagonal")
rv1

From the above result, we can see that most of the cell-type-specific TFs only at activity level fall into fetal liver HSCs and adult bone marrow HSCs. Next, to further explore the disribution of these cell-type-specific TFs, we colored the points with their cell-type-specificity label.

rv2 <- RadvizEnrichPlot(x = cts_act_frm, color_by='stage', anchors = stages, optim_anchor = F, recenter = stages, text_size = 4)
rv2

In the above figure, we found that the cell-type-specificity label showed relative high correlation with their location. Note that the rv generate from RadvizEnrichPlot() function is a ggplot object, which means that you can modify this object with functions in ggplot2 package (e.g. adding gene labels).

candidates <- c("Sox6","Ddit3","Cebpb")
names(candidates) <- candidates
cts_act_frm$gene_label <- candidates[cts_act_frm$regulon]

rv2 <- RadvizEnrichPlot(x = cts_act_frm, color_by='stage', anchors = stages, optim_anchor = F, recenter = stages, text_size = 4)+
  geom_text(mapping = aes(label=gene_label),na.rm = TRUE)
rv2

Besides, we noticed that some of the above TFs showed cell-type-specificity in more than one cell types. To better describe cell-type-specific distribution among these cell types, we can also use Venn diagram to visualize these TFs:

library(UpSetR)

## 
## 载入程辑包：'UpSetR'

## The following object is masked from 'package:lattice':
## 
##     histogram

cts_act_only_stat <- table(cts_act_only$regulon,cts_act_only$stage)
attributes(cts_act_only_stat)$class <- "matrix"
data_set <- as.data.frame(cts_act_only_stat)
data_set <- data_set[,stages]

setR <- upset(data_set, nsets = 6, nintersects = 40, mb.ratio = c(0.5, 0.5),
      order.by = c("degree","freq"), decreasing = c(F,TRUE), point.size=5,
      line.size=2)
setR